      program main_htran
	implicit none
!-------------------------------------------------------------------------
! Transform geodetic height to orthometric/normal height, or vice-visa,
! the geoid height or height anomaly are interpolated from a 1 x 1 minute 
! geoid/quasi-geoid model
!-------------------------------------------------------------------------
      character*200    fmod,fin,fout
      integer*4        imod,iin,iout
      integer*4        nrow,ncol,i,j,iwind,ios
      integer*4        itrs,ifmt,prec,ierr,utm(3)
	double precision minlat,maxlat,minlon,maxlon,lons,lats
      double precision dlat,dlon,und
	real*4,allocatable:: grid(:,:)
	double precision,allocatable:: codt(:,:),h1(:),h2(:)
      character(len=200) :: arg,utm_str
      integer*4          :: num_args,pos1,pos2
      real*4::Ellipsoid_a,Ellipsoid_f
      real*8 ::X, Y,lat,lon,CENTRAL_LON

      INTEGER LINE_NUMBER,GETFILE_LINENUMBER
!-------------------------------------------------------------------------
        ! 获取命令行参数
      num_args = command_argument_count()
      if (num_args < 5) then
            print *, "Example:"
            stop
      end if

      call get_command_argument(1, fmod)  ! 模型文件路径
      call get_command_argument(2, fin)   ! 输入文件路径
      call get_command_argument(3, fout)  ! 输出文件路径

       ! 解析转换类型 itrs
      call get_command_argument(4, arg) !-1是大地高-正高  1是正高-大地高
      read(arg, *, iostat=ierr) itrs
      if (ierr /= 0 .or. (itrs /= -1 .and. itrs /= 1)) then
        print *, "Error: itrs must be -1 or 1."
        stop
      end if

      ! 解析UTM参数
      call get_command_argument(5, utm_str) !(0，0，0 )  
      pos1 = index(utm_str, ',')
      pos2 = index(utm_str(pos1+1:), ',') + pos1
      if (pos1 == 0 .or. pos2 == pos1) then
        print *, "Error: Invalid UTM format. Expected 'utm1,utm2,utm3'"
        stop
      end if
      read(utm_str(1:pos1-1), *, iostat=ierr) utm(1) !第一个参数代表坐标系  0是大地坐标系  1是utm -1是高斯
      read(utm_str(pos1+1:pos2-1), *, iostat=ierr) utm(2) !如果是1   后边两个分别代表半球和带号
      read(utm_str(pos2+1:), *, iostat=ierr) utm(3) !带号
      if (ierr /= 0) then
        print *, "Error: UTM parameters must be integers."
        stop
      end if
      imod=31
	iin=32
      iout=51
	print *
      print *,'========== Start of height transformation =========='
c c
c
c Geoid model file
c
c	fmod='angola_quasi_geoid.grd'
c
c Input height file
c
c	fin='SC33E021022DEMU_geoh.blh'
c     fin='SC33E021022DEMU_geoh.xyh'
c
c Onput height file
c
c      fout='SC33E021022DEMU_orth_out.blh'
c      fout='SC33E021022DEMU_orth_out.xyh'
c
c Type of height transformation
c
c      itrs=-1 ! geodetic height           ==> normal/orthometric height
c	itrs=1  ! normal/orthometric height ==> geodetic height
	if(itrs.ne.-1.and.itrs.ne.1) then
	ierr=1
	call reperr(ierr)
	stop
	endif
c
c Type of coordinate
c  utm(3): geodetic coordinate    if utm(1)=0
c          UTM coordinate         if utm(1)=1
c          utm(2): Hemisphere option, 1 for north, -1 for south
c          utm(3): UTM Zone number
c          utm(2), utm(3) only working if utm(1)=1
c	utm(1)=1
c	utm(2)=-1
c	utm(3)=33
	if(utm(1).ne.0.and.utm(1).ne.1.and.utm(1).ne.-1) then
	ierr=2
	call reperr(ierr)
	stop
	endif
	if(utm(1).eq.1) then
	if((utm(2).ne.1.and.utm(2).ne.-1).or.
     $   (utm(3).lt.1.or.utm(3).gt.61)) then
	ierr=3
c	print *,utm(1),utm(2),utm(3)
	call reperr(ierr)
	stop
	endif
	endif
c
c Format of input and output height file
c
      ifmt=0 ! ASCII format
c	ifmt=1 ! IMG format
c	ifmt=2 ! GeoTiff format

c
c Spline window size for interpolation
c     
      iwind=4
!
!     Read geoid model in quasi SRTM binary grid
!
      open(imod,file=fmod,status='old',err=601,form='unformatted',
     $     iostat=ios)
	if(ios.gt.0) goto 601 ! ierr=4
	
	if(ifmt.eq.0) then

      open(iin,file=fin,status='old',err=602,iostat=ios)
	if(ios.gt.0) goto 602 ! ierr=5

	open(iout,file=fout,status='replace',err=603,iostat=ios)
	if(ios.gt.0) goto 603 ! ierr=6

	endif
	
      read(imod,iostat=ios,err=604) minlat,maxlat,minlon,maxlon,lats,  ! ierr=7
     $                              lons,nrow,ncol
c      print *,'binary grid file header: ',
c     $        minlat,maxlat,minlon,maxlon,lats,lons,nrows,ncols
      if(minlon.gt.180.d0) minlon=minlon-360
      if(maxlon.gt.180.d0) maxlon=maxlon-360
      allocate(grid(nrow,ncol))
      read(imod,iostat=ios,err=604) ((grid(i,j),j=1,ncol),i=nrow,1,-1)  ! ierr=7
      close(imod)

	i=0
14    continue
      read(iin,*,iostat=ios,err=605,end=15) 
	if(ios.gt.0) goto 605 ! ierr=8
      i=i+1
      goto 14
15    continue
      prec=i  
      allocate(codt(prec,2))
	allocate(h1(prec))
	allocate(h2(prec))
    	
	rewind(iin)
c! Start of loop i
      do i=1,prec
      read(iin,*,iostat=ios,err=605) codt(i,1),codt(i,2),h1(i)
	if(ios.gt.0) goto 605  ! ierr=8
	if(utm(1).eq.0) then
	if(dabs(codt(i,1)).gt.90.d0.or.
     $   codt(i,2).lt.-180.d0.or.codt(i,2).gt.360.d0) then
	ierr=9
	call reperr(ierr)
	stop
	endif
	endif
      if(utm(1).eq.0) then
	dlat=codt(i,1)
	dlon=codt(i,2)
	endif
c Transofrm UTM coordinates to geodetic coordinates, if utm(1)=1
      if(utm(1).eq.1) then 
      call ugtran(0,utm(2),utm(3),codt(i,1),codt(i,2),
     $                dlat,dlon)
	endif
      if(dlon.gt.180.d0) then 
      dlon=dlon-360
      call interp(iwind,0.d0,grid,minlat,minlon,lats,lons,nrow,ncol,
     $            nrow,ncol,dlat,dlon,und)
      h2(i)=h1(i)+itrs*und
      write(iout,300) codt(i,1),codt(i,2),h2(i)
	endif

ccTransofrm GUss coordinates to geodetic coordinates, if utm(1)=-1
      if(utm(1).eq.-1) then 
      Ellipsoid_a=6378245.0d0
      Ellipsoid_f=1.0d0 / 298.3d0
      CENTRAL_LON=utm(3)*6-3
c X=3840570.748385d0
c Y=749372.187674d0
      X=codt(i,1)
      Y=codt(i,2)
      call gauss_inverse(X, Y, CENTRAL_LON, 6378245.0d0, 
     $ 1.0d0 / 298.3d0, lat, lon)
      print *, "反算结果：纬度=", lat, ", 经度=",lon
      write(iout,300) lat,lon,h1(i)
      endif

      enddo ! End of loop i
      close(iin)

c do i=1,prec
c write(iout,300) codt(i,1),codt(i,2),h2(i)
c enddo
	close(iout)

      deallocate(grid)
	deallocate(codt)
	deallocate(h1)
	deallocate(h2)
	print *
	print *,'========== End of height transformation =========='
	print *

300   format(2f30.20,f20.10)

      goto 701

601   ierr=4
      goto 700
602   ierr=5
      goto 700
603   ierr=6
      goto 700
604   ierr=7
      goto 700
605   ierr=8
      goto 700 

700   continue
      call reperr(ierr)
701   continue

	stop
	end

      subroutine ugtran(ityp,hemi,zone,easting,northing,dlat,dlon)
	implicit none
c  Input variables
      integer*4        ityp,hemi,zone
	double precision easting,northing ! Input if ityp=0; Output if ityp=1
c  Output variables
      double precision dlat,dlon ! Input if ityp=1; Output if ityp=0
c  Local variables
        double precision pi,dtr
        double precision a,f,e2,e4,e6,k0,ka,easting0,northing0
        double precision tmp1,tmp2,tmp3,tmp4
        double precision rlat,rlon,dlon0,rlon0,east,north
        double precision v,s,t,c,biga,biga2,biga3,biga4,biga5,biga6
        double precision sinlat,coslat,tanlat,sinf,cosf,tanf
        double precision e1,e11,e12,e13,e14,ee2,j1,j2,j3,j4,c1,c12
        double precision m,mu,fp,buf1,r1,n1,d,d2,d3,d4,d5,d6,t1,t12
        double precision q1,q2,q3,q4,q6,q7

        pi=4.d0*datan(1.d0)
        dtr=pi/180.d0
c
c  Parameters of WGS84 geodetic reference system
c
      a=6378137.d0
      f=1.d0/298.257223563d0 ! tide free
      e2=0.00669437999014d0
      e4=e2*e2
      e6=e4*e2

      k0=0.9996d0
      ka=k0*a
      easting0=500000.d0

      if(ityp.eq.0) then ! Start of ityp=0

      tmp1=a*(1.d0-e2/4.d0-3.d0*e4/64.d0-5.d0*e6/256.d0)
      tmp2=dsqrt(1.d0-e2)
      tmp3=a*(1.d0-e2)
      e1=(1.d0-tmp2)/(1.d0+tmp2)
      e12=e1*e1
      e13=e12*e1
      e14=e13*e1
      j1=(1.5d0*e1-27.d0*e13/32.d0)
      j2=(21.d0*e12/16.d0-55.d0*e14/32.d0)
      j3=151.d0*e13/96.d0
      j4=1097.d0*e14/512.d0
      ee2=e2/(1.d0-e2)

c     if(easting.lt.167000.d0.or.easting.gt.833000.d0) then
c      print *,'Wrong easting, stop !'
c      stop
c     endif
c      if(northing.lt.0.d0) then
c      print *,easting,northing,h
c      northing=10000000.d0+northing
c      print *,easting,northing,h
c      print *,'Wrong northing, stop !'
c      stop
c      endif
      if(hemi.eq.1)  northing0=0.d0
      if(hemi.eq.-1) northing0=10000000.d0
      north=northing-northing0
      east=easting-easting0
      m=north/k0
      mu=m/tmp1
      fp=mu+j1*dsin(2.d0*mu)+j2*dsin(4.d0*mu)+j3*dsin(6.d0*mu)+
     $      j4*dsin(8.d0*mu)
       sinf=dsin(fp)
       cosf=dcos(fp)
       tanf=dtan(fp)
       c1=ee2*cosf*cosf
       c12=c1*c1
       t1=tanf*tanf
       t12=t1*t1
        buf1=dsqrt(1.d0-e2*sinf*sinf)
        r1=tmp3/buf1/buf1/buf1
        n1=a/buf1
        d=east/n1/k0
        d2=d*d
        d3=d2*d
        d4=d3*d
        d5=d4*d
        d6=d5*d
        q1=n1*tanf/r1
        q2=d2/2.d0
        q3=(5.d0+3*t1+10.d0*c1-4.d0*c12-9.d0*ee2)*d4/24.d0
        q4=(61.d0+90.d0*t1+298.d0*c1+45.d0*t12-3.d0*c12-252.d0*ee2)*
     $    d6/720.d0
        q6=(1.d0+2.d0*t1+c1)*d3/6.d0
        q7=(5.d0-2.d0*c1+28.d0*t1-3.d0*c12+8.d0*ee2+24.d0*t12)*d5/120.d0
        dlon0=(zone-1)*6-180+3
        dlat=(fp-q1*(q2-q3+q4))/dtr
        dlon=dlon0+(d-q6+q7)/cosf/dtr

      endif ! End of ityp=0


      if(ityp.eq.1) then ! Start of ityp=1

      tmp1=1.d0-e2/4.d0-3.d0*e4/64.d0-5.d0*e6/256.d0
      tmp2=3.d0*e2/8.d0+3.d0*e4/32.d0+45.d0*e6/1024.d0
      tmp3=15.d0*e4/256.d0+45.d0*e6/1024.d0
      tmp4=35.d0*e6/3072.d0

      if(dlat.ge.0.d0) then
       northing0=0.d0
       hemi=1
      else
       northing0=10000000.d0
        hemi=-1
      endif
        if(dlon.gt.180.d0) dlon=dlon-360.d0 
        if(dlon.lt.0.d0) then
        zone=aint(dlon/6.d0)+30
        else
        zone=aint(dlon/6.d0)+31
        endif
        rlat=dlat*dtr
        rlon=dlon*dtr
        sinlat=dsin(rlat)
        coslat=dcos(rlat)
        tanlat=dtan(rlat)
        v=1.d0/dsqrt(1.d0-e2*sinlat*sinlat)
        rlon0=((zone-1)*6-180+3)*dtr
c	print *,(zone-1)*6-180+3
        biga=(rlon-rlon0)*coslat
        biga2=biga*biga
        biga3=biga2*biga
        biga4=biga3*biga
        biga5=biga4*biga
        biga6=biga5*biga
        s=tmp1*rlat-tmp2*dsin(2.d0*rlat)+tmp3*dsin(4.d0*rlat)-
     $  tmp4*dsin(6.d0*rlat)
        t=tanlat*tanlat
        c=e2*coslat*coslat/(1-e2)
        easting=easting0+ka*v*(biga+(1-t+c)*biga3/6.d0+
     $        (5.d0-18.d0*t+t*t)*biga5/120.d0)
      northing=northing0+ka*(s+v*tanlat*(biga2/2.d0+
     $         (5.d0-t+9.d0*c+4.d0*c*c)*biga4/24.d0+
     $         (61.d0-58.d0*t+t*t)*biga6/720.d0))
c      write(iout,300) hemi,zone,northing,easting

      endif ! End of ityp=1

      return
      end

!
!-----------------------------------------------------------------------
!
      SUBROUTINE INTERP(iwO,DMIN,H,PHIS,DLAW,DDFI,DDLA,NPHI,NDLA,
     $                  IPDIM,ILDIM,PHI,DLA,VALINT)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!                                                                      C
!     SUBROUTINE FOR INTERPOLATION OF VALUES FROM A STANDARD DTM-GRID  C
!     TO INDIVIDUAL STATION LOCATIONS.                                 C
!                                                                      C
!                                                                      C
!     INPUT PARAMETERS...                                              C
!     ===================                                              C
!     iwO...    A SPLINE WINDOW OF SIZE 'iwO' X 'iwO' WILL BE C
!                  USED AROUND EACH STATION. IF 'iwO' IS 0 OR 1,    C
!                  BILINEAR INTERPOLATION WILL BE USED.                C
!     DMIN...      MINIMUM ACCEPTABLE DISTANCE FROM THE GRID EDGE IN   C
!                  KM (USEFUL FOR FFT GRIDS).                          C
!     H...         2D DATA ARRAY (ELEMENT (1,1) IN SW CORNER).         C
!     PHIS,DLAW... LATITUDE AND LONGITUDE OF SW GRID POINT.            C
!     DDFI,DDLA... GRID SPACING IN LATITUDE AND LONGITUDE DIRECTION.   C
!     NPHI,NDLA... NUMBER OF GRID POINTS IN LATITUDE AND LONGITUDE     C
!                  DIRECTION.                                          C
!     IPDIM,ILDIM..DIMENSIONS OF 2D DATA ARRAY 'H' AS DECLARED IN THE  C
!                  CALLING PROGRAM.                                    C
!     PHI,DLA...   LATITUDE AND LONGITUDE OF INTERPOLATION POINT.      C
!                                                                      C
!                                                                      C
!     OUTPUT PARAMETERS...                                             C
!     ====================                                             C
!     VALINT...    INTERPOLATED VALUE.                                 C
!                                                                      C
!                                                                      C
!     EXECUTION TIME ON CDC 990 IS...                                  C
!     ===============================                                  C
!     +------------------+-------------------+-------------------+     C
!     I  INTERPOLATION   I  OPT=LOW          I  OPT=HIGH         I     C
!     I------------------I-------------------I-------------------I     C
!     I  BILINEAR        I  1.44 MSEC/STAT.  I  1.44 MSEC/STAT.  I     C
!     I  3 X 3 SPLINE    I  1.53 MSEC/STAT.  I  1.51 MSEC/STAT.  I     C
!     I  5 X 5 SPLINE    I  1.70 MSEC/STAT.  I  1.67 MSEC/STAT.  I     C
!     I  7 X 7 SPLINE    I  2.02 MSEC/STAT.  I  1.74 MSEC/STAT.  I     C
!     I  9 X 9 SPLINE    I  2.31 MSEC/STAT.  I  2.00 MSEC/STAT.  I     C
!     +------------------+-------------------+-------------------+     C
!                                                                      C
!                                                                      C
!     PROGRAM CREATION BY...   H. DENKER          MAY 30, 1987         C
!                              H. DENKER          MARCH 13, 1989       C
!                                                                      C
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      PARAMETER(IPA1=20)
      IMPLICIT REAL*8(A-H,O-Z)
      LOGICAL LODD
      REAL*4 H(IPDIM,ILDIM)
c      REAL*8 H(IPDIM,ILDIM)
      DIMENSION A(IPA1),R(IPA1),Q(IPA1),HC(IPA1)
      IDIM1=IPA1
       TWOPI=6.28318530717959D0
         RHO=360.D0/TWOPI
      REARTH=6371000.D0
      IF(iwO.LT.2) iwO=2
      IF(iwO.GT.IDIM1) iwO=IDIM1
      ILIM=DMIN*1000.*RHO/(REARTH*DDFI)
      JLIM=DMIN*1000.*RHO/(REARTH*DDLA*COS((PHIS+DDFI*NPHI/2.)/RHO))
      LODD=(iwO/2)*2.NE.iwO
      RI=(PHI-PHIS)/DDFI
      RJ=(DLA-DLAW)/DDLA
      IF(LODD) THEN
        I0=RI-0.5
        J0=RJ-0.5
      ELSE
        I0=RI
        J0=RJ
      ENDIF
      I0=I0-iwO/2+1
      J0=J0-iwO/2+1
      II=I0+iwO-1
      JJ=J0+iwO-1
      IF(I0.LT.0 .OR. II.GE.NPHI .OR. J0.LT.0 .OR. JJ.GE.NDLA) THEN
        WRITE(6,7008) PHI,DLA
        VALINT=999999.
        RETURN
      ELSEIF(I0.LT.ILIM .OR. II.GT.NPHI-ILIM .OR. J0.LT.JLIM .OR.
     $  JJ.GT.NDLA-JLIM) THEN
        IF(NPOINT.LE.ILIST) WRITE(6,7009) PHI,DLA
        VALINT=999999.
        RETURN
      ENDIF
7008  FORMAT(' ',2F12.6,' STATION TOO NEAR GRID BOUNDARY  - NO INT.'
     .,' POSSIBLE|')
7009  FORMAT(' ',2F12.6,' STATION OUTSIDE ACCEPTABLE AREA - NO INT.'
     .,' PERFORMED|')
      IF(iwO.GT.2) THEN
        DO 110 I=1,iwO
          DO 111 J=1,iwO
            A(J)=H(I0+I,J0+J)
111       CONTINUE
          CALL INITSP(A,iwO,R,Q)
          HC(I)=SPLINE(RJ-J0+1.,A,iwO,R)
110     CONTINUE
        CALL INITSP(HC,iwO,R,Q)
        VALINT=SPLINE(RI-I0+1.,HC,iwO,R)
      ELSE
        VALINT=BILIN(RI+1.,RJ+1.,H,NPHI,NDLA,IPDIM,ILDIM)
      ENDIF
      RETURN
      END
!
!-----------------------------------------------------------------------
!
      FUNCTION BILIN(RI,RJ,A,IMAX,JMAX,IADIM,JADIM)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!                                                                      C
!                           B I L I N                                  C
!                                                                      C
!  INTERPOLATES VALUES IN AN ARRAY A USING BILINEAR                    C
!  (PARABOLIC HYPERBOLOID) INTERPOLATION.                              C
!                                                                      C
!----------------------------------------------------------------------C
!                                                                      C
!  PARAMETERS:                                                         C
!                                                                      C
!  BILIN...       INTERPOLATED VALUE                                   C
!                                                                      C
!  RI, RJ...      INTERPOLATION ARGUMENT, (1,1) IN LOWER LEFT CORNER,  C
!                 (IMAX, JMAX) IN UPPER RIGHT.                         C
!                                                                      C
!  A...           INTEGER*2 ARRAY WITH ARGUMENTS                       C
!                                                                      C
!  IMAX, JMAX...  NUMBER OF POINTS IN GRID                             C
!                                                                      C
!  IADIM, JADIM...DECLARED DIMENSIONS OF 'A'                           C
!                                                                      C
!  OUTSIDE AREA COVERED BY 'A' THE FUNCTION RETURNS THE VALUE OF       C
!  THE NEAREST BOUNDARY POINT.                                         C
!                                                                      C
!----------------------------------------------------------------------C
!                                                                      C
!  PROGRAMMER:                                                         C
!  RENE FORSBERG, JULY 1983                                            C
!                                                                      C
!  MODIFICATIONS BY:                                                   C
!  HEINER DENKER, 07/01/1987                                           C
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IMPLICIT REAL*8(A-H,O-Z)
c      DIMENSION A(IADIM, JADIM)
      REAL*4 A(IADIM, JADIM)
      IN = IFRAC(RI)
      IE = IFRAC(RJ)
      RN = RI - IN
      RE = RJ - IE
      IF (IN.LT.1) THEN
        IN = 1
        RN = 0.0
      ELSEIF (IN.GE.IMAX) THEN
        IN = IMAX-1
        RN = 1.0
      ENDIF
      IF (IE.LT.1) THEN
        IE = 1
        RE = 0.0
      ELSEIF (IE.GE.JMAX) THEN
        IE = JMAX-1
        RE = 1.0
      ENDIF
      RNM1=1.-RN
      REM1=1.-RE
      BILIN = RNM1*REM1*A(IN,IE) +
     .RN*REM1*A(IN+1,IE) + RNM1*RE*A(IN,IE+1) +
     .RN*RE*A(IN+1,IE+1)
      RETURN
      END
!
!-----------------------------------------------------------------------
!
      SUBROUTINE INITSP(Y, N, R, Q)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!                                                                      C
!                      I N I T S P                                     C
!                                                                      C
!  INITIALIZATION PROCEDURE FOR FAST 1-DIMENSIONAL EQUIDISTANT         C
!  SPLINE INTERPOLATION, WITH FREE BOUNDARY END CONDITIONS             C
!  REFERENCE: JOSEF STOER: EINFUHRUNG IN DIE NUMERISCHE MATHEMATIK     C
!  I, SPRINGER 1972, PAGE 82 AND 86.                                   C
!                                                                      C
!----------------------------------------------------------------------C
!                                                                      C
!  PARAMETERS (REAL):                                                  C
!                                                                      C
!  Y...   GIVEN VALUES, Y(1), ..., Y(N)                                C
!                                                                      C
!  R...   SPLINE MOMENTS (1 ... N), TO BE USED BY FUNCTION 'SPLINE'    C
!                                                                      C
!  Q...   WORK-ARRAY, DECLARED AT LEAST 1:N                            C
!                                                                      C
!----------------------------------------------------------------------C
!                                                                      C
!  RENE FORSBERG, JULY 1983                                            C
!                                                                      C
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION Y(N), R(N), Q(N)
      Q(1) = 0.0
      R(1) = 0.0
      DO 11 K = 2, N-1
        P = Q(K-1)/2+2
        Q(K) = -0.5/P
        R(K) = (3*(Y(K+1)-2*Y(K)+Y(K-1)) - R(K-1)/2)/P
   11 CONTINUE
      R(N) = 0.0
      DO 12 K = N-1, 2, -1
        R(K) = Q(K)*R(K+1)+R(K)
   12 CONTINUE
      RETURN
      END
!
!-----------------------------------------------------------------------
!
      FUNCTION SPLINE(X, Y, N, R)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!                                                                      C
!                          S P L I N E                                 C
!                                                                      C
!  FAST ONE-DIMENSIONAL EQUIDISTANT SPLINE INTERPOLATION FUNCTION.     C
!  REFERENCE: JOSEF STOER: EINFUHRUNG IN DIE NUMERISCHE MATHEMATIK     C
!  I, SPRINGER 1972, PAGE 81.                                          C
!                                                                      C
!----------------------------------------------------------------------C
!                                                                      C
!  PARAMETERS:                                                         C
!                                                                      C
!  X...  INTERPOLATION ARGUMENT (REAL), X = 1 FIRST DATA-POINT,        C
!        X = N LAST DATA-POINT. OUTSIDE THE RANGE LINEAR EXTRA-        C
!        POLATION IS USED.                                             C
!                                                                      C
!  Y...  REAL*8 ARRAY, 1 .. N : DATA VALUES                            C
!                                                                      C
!  R...  DO: SPLINE MOMENTS CALCULATED BY SUBROUTINE 'INITSP'          C
!                                                                      C
!----------------------------------------------------------------------C
!                                                                      C
!  PROGRAMMER:                                                         C
!  RENE FORSBERG, JUNE 1983                                            C
!                                                                      C
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IMPLICIT REAL*8(A-H, P-Z)
      DIMENSION Y(N), R(N)
      IF (X.LT.1) THEN
        SPLINE = Y(1) + (X-1)*(Y(2)-Y(1)-R(2)/6)
      ELSEIF (X.GT.N) THEN
        SPLINE = Y(N) + (X-N)*(Y(N)-Y(N-1)+R(N-1)/6)
      ELSE
        J = IFRAC(X)
        XX = X - J
        SPLINE = Y(J) +
     .           XX * ((Y(J+1)-Y(J)-R(J)/3-R(J+1)/6) +
     .           XX * (R(J)/2 +
     .           XX * (R(J+1)-R(J))/6))
      ENDIF
      RETURN
      END
!
!-----------------------------------------------------------------------
!
      FUNCTION IFRAC(R)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!                                                                      C
!                    F U N C T I O N   I F R A C                       C
!                  ===============================                     C
!                                                                      C
!  SUBROUTINE GIVING TRUE INTEGER PART OF A REAL E.G.                  C
!                                                                      C
!    FOR   1. = R < 2.   IS    IFRAC = 1                               C
!    FOR   0. = R < 1.   IS    IFRAC = 0                               C
!    FOR  -1. = R < 0.   IS    IFRAC =-1                               C
!    FOR  -2. = R <-1.   IS    IFRAC =-2                               C
!                                                                      C
!  RF, JUNE 1983                                                       C
!  HD, JANUARY 1987                                                    C
!                                                                      C
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IMPLICIT REAL*8(A-H,O-Z)
      IFRAC=R
      IF (R.GE.0) RETURN
      IF (R.EQ.IFRAC) RETURN
      IFRAC = IFRAC - 1
      RETURN
      END
!
!-----------------------------------------------------------------------
!

	subroutine reperr(ierr)
	implicit none

	integer*4 ierr

	print *
	print *,'========== Report error in htran =========='
	print *

	if(ierr.eq.1)
     $write(*,301) 'ierr = ',ierr,'(Height transformation type)'
	if(ierr.eq.2)
     $write(*,301) 'ierr = ',ierr,'(Coordinate type)'
	if(ierr.eq.3)
     $write(*,301) 'ierr = ',ierr,'(UTM hemisphere and zone number)'
	if(ierr.eq.4)
     $write(*,301) 'ierr = ',ierr,'(Geoid model file opening)'
	if(ierr.eq.5)
     $write(*,301) 'ierr = ',ierr,'(Input height file opening)'
	if(ierr.eq.6)
     $write(*,301) 'ierr = ',ierr,'(Output height file opening)'
	if(ierr.eq.7)
     $write(*,301) 'ierr = ',ierr,'(Geoid model file reading)'
	if(ierr.eq.8)
     $write(*,301) 'ierr = ',ierr,'(Input height file reading)'
	if(ierr.eq.9)
     $write(*,301) 'ierr = ',ierr,'(Input coordinate)'
	
	print *
	print *

301   format(a8,i2,a35)

	return
	end

      SUBROUTINE GAUSS_INVERSE(X, Y, CENTRAL_LON, A, F, LAT, LON)
C 显式声明变量
      IMPLICIT NONE
      real*8 X, Y, CENTRAL_LON, A, F, LAT, LON
      REAL*8 Y_, L0, E2, E4, E6, K0, K2, K4, K6, BF, SIN_BF, COS_BF
      REAL*8 TAN_BF, NF, ETA2, RHO, T, X_, Y2, Y4, B, L
      INTEGER I
C 去除 500 公里偏移
      Y_ = Y - 500000.0D0
      L0 = CENTRAL_LON * ATAN(1.0D0) / 45.0D0
C 这里直接传入 a 和 f，
      E2 = F * 2 - F * F
      E4 = E2**2
      E6 = E2**3
C 计算 K0, K2, K4, K6
      K0 = 1.0D0 + (3.0D0 / 4.0D0) * E2 + (45.0D0 / 64.0D0) * 
     $      E2**2 + (175.0D0 / 256.0D0) * E2**3
      K2 = -(3.0D0 / 8.0D0) * E2 - (15.0D0 / 32.0D0) * E2**2 - 
     $      (525.0D0 / 1024.0D0) * E2**3
      K4 = (15.0D0 / 256.0D0) * E2**2 + (105.0D0 / 1024.0D0) * E2**3
      K6 = -(35.0D0 / 3072.0D0) * E2**3
C 计算底点纬度 Bf（迭代法）
      BF = X / (A * (1.0D0 - E2 / 4.0D0 - 3.0D0 * E4 / 64.0D0 - 
     $      5.0D0 * E6 / 256.0D0))
      DO 10 I = 1, 100
         SIN_BF = SIN(BF)
         COS_BF = COS(BF)
         TAN_BF = TAN(BF)
         NF = A / SQRT(1.0D0 - E2 * SIN_BF**2)
         ETA2 = E2 * COS_BF**2 / (1.0D0 - E2)
         RHO = A * (1.0D0 - E2) / (1.0D0 - E2 * SIN_BF**2)**(1.5D0)
         BF = BF + (X - (A * (1.0D0 - E2) * (K0 * BF + K2 * 
     $      SIN(2.0D0 * BF) + K4 * SIN(4.0D0 * BF) + K6 * 
     $      SIN(6.0D0 * BF)))) / RHO
         IF (ABS((X - (A * (1.0D0 - E2) * (K0 * BF + K2 * 
     $      SIN(2.0D0 * BF) + K4 * SIN(4.0D0 * BF) + K6 * 
     $      SIN(6.0D0 * BF)))) 
     $      / RHO) < 1.0D-12) GOTO 20
 10   CONTINUE
 20   CONTINUE
C 计算辅助量
      SIN_BF = SIN(BF)
      COS_BF = COS(BF)
      TAN_BF = TAN(BF)
      NF = A / SQRT(1.0D0 - E2 * SIN_BF**2)
      ETA2 = E2 * COS_BF**2 / (1.0D0 - E2)
      RHO = A * (1.0D0 - E2) / (1.0D0 - E2 * SIN_BF**2)**(1.5D0)
      T = TAN_BF
      X_ = X
      Y2 = Y_**2
      Y4 = Y_**4
C 计算 B 和 l
      B = BF - (Y2 * T) / (2.0D0 * RHO * NF) * (1.0D0 - (Y2) / 
     $      (12.0D0 * NF**2) * (5.0D0 + 3.0D0 * T**2 + ETA2 - 9.0D0
     $       * ETA2 * T**2))
      L = (Y_) / (NF * COS_BF) * (1.0D0 - (Y2) / (6.0D0 * NF**2) * 
     $      (1.0D0 + 2.0D0 * T**2 + ETA2) + (Y4) / (120.0D0 * NF**4) 
     $      * (5.0D0 + 28.0D0 * T + 24.0D0 * T**2 + 6.0D0 * ETA2 + 8.0D0 
     $      * T * ETA2))
C 转换为角度
      LAT = B * 45.0D0 / ATAN(1.0D0)
      LON = L * 45.0D0 / ATAN(1.0D0) + CENTRAL_LON
      RETURN
      END

c       INTEGER FUNCTION GETFILE_LINENUMBER(FILE_NAME, IFILEUNIT)
c C 显式声明变量
c       CHARACTER*(*) FILE_NAME
c       INTEGER IFILEUNIT
c       INTEGER IOS, IO_STATUS
c       CHARACTER*1 CDUMMY
c       CHARACTER*100 ERROR_MESSAGE
c C 初始化返回值
c       GETFILE_LINENUMBER = 0
c C 打开文件
c       OPEN(UNIT = IFILEUNIT, FILE = FILE_NAME, STATUS = 'OLD',
c      &       ACTION = 'READ', IOSTAT = IO_STATUS)
c       IF (IO_STATUS.NE. 0) THEN
c          WRITE(*,*) 'Error opening file: ', ERROR_MESSAGE
c          STOP
c       END IF
c c C 将文件指针重置到文件开头
c c       REWIND(IFILEUNIT)
c c C 循环读取文件行
c  10   READ(IFILEUNIT, *, IOSTAT = IOS) CDUMMY
c       IF (IOS.EQ. 0) THEN
c          GETFILE_LINENUMBER = GETFILE_LINENUMBER + 1
c          GOTO 10
c       END IF
c       RETURN
c       END FUNCTION GETFILE_LINENUMBER